﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "pqeuler2")]
    public static unsafe void pqeuler2(double t, double h, IntPtr y_ptr, IntPtr z_ptr, double eps, IntPtr f_abc_ptr)
    {
        double* y = (double*)y_ptr.ToPointer();
        double* z = (double*)z_ptr.ToPointer();
        f_abc = Marshal.GetDelegateForFunctionPointer<delegatefunc_abc>(f_abc_ptr);

        pqeuler2(t, h, y, z, eps);
    }

    /// <summary>
    /// 连分式法求解二阶初值
    /// f计算二阶微分方程的右端函数f(t,y,z)。
    /// </summary>
    /// <param name="t">自变量起点值</param>
    /// <param name="h">步长</param>
    /// <param name="y">存放函数初值。返回终点函数值。</param>
    /// <param name="z">存放函数一阶导数初值。返回终点函数一阶导数值。</param>
    /// <param name="eps">精度要求</param>
    public static unsafe void pqeuler2(double t, double h, double* y, double* z, double eps)
    {
        int j, il, flag, m;
        double yy, zz;
        double h0, d, ys0, ys1 = 0.0, zs = 0.0, y0, z0;
        double* yb = stackalloc double[10];
        double* zb = stackalloc double[10];
        double* hh = stackalloc double[10];
        double* gy = stackalloc double[10];
        double* gz = stackalloc double[10];

        yy = *y;
        zz = *z;
        il = 0;
        flag = 0;
        m = 1;
        h0 = h;

        // Euler方法计算初值gy[0]与gz[0]
        euler21(t, h0, &yy, &zz, m);
        y0 = yy;
        z0 = zz;
        while ((il < 20) && (flag == 0))
        {
            il = il + 1;
            hh[0] = h0;
            gy[0] = y0;
            gz[0] = z0;
            // 计算yb[0]与zb[0]
            yb[0] = gy[0];
            zb[0] = gz[0];
            j = 1;
            ys1 = gy[0];
            while (j <= 7)
            {
                yy = *y;
                zz = *z;
                m = m + m;
                hh[j] = hh[j - 1] / 2.0;
                //Euler方法计算新近似值gy[j]与gz[j]
                euler21(t, hh[j], &yy, &zz, m);
                gy[j] = yy;
                gz[j] = zz;
                //计算yb[j]
                funpqj(hh, gy, yb, j);
                //计算zb[j]
                funpqj(hh, gz, zb, j);
                ys0 = ys1;
                //连分式法计算积分近似值ys1
                ys1 = funpqv(hh, yb, j, 0.0);
                //连分式法计算积分近似值zs
                zs = funpqv(hh, zb, j, 0.0);
                d = Math.Abs(ys1 - ys0);
                if (d >= eps) j = j + 1;
                else j = 10;
            }
            h0 = hh[j - 1];
            y0 = gy[j - 1];
            z0 = gz[j - 1];
            if (j == 10)
            {
                flag = 1;
            }
        }
        *y = ys1;
        *z = zs;
        return;
    }

    /*
    // 求解二阶初值连分式法例
      int main()
      {
          int j;
          double t,h,eps,y,z;
          double  pqeuler2f(double, double, double);
          y=0.0; z=0.701836;
          t=0.0; h=0.1; eps=0.0000001;
          cout <<"t = " <<setw(6) <<t;
          cout <<setw(6) <<"y = " <<setw(10) <<y;
          cout <<setw(6) <<"z = " <<setw(10) <<z;
          cout <<endl;
          for (j=1; j<=10; j++)
          { 
              pqeuler2(t,h,&y,&z,eps,pqeuler2f);
              t=t+h;
              cout <<"t = " <<setw(6) <<t;
              cout <<setw(6) <<"y = " <<setw(10) <<y;
              cout <<setw(6) <<"z = " <<setw(10) <<z;
              cout <<endl;
          }	  
          return 0;
      }
    // 计算二阶微分方程的右端函数f(t,y,z)
      double pqeuler2f(double t, double y, double z)
      { 
          double d;
          d=t+y;
          return(d);
      }
    */
}

